/*
 * Example_CO2Solubility.h
 *
 *  Created on: 12 Sep 2011
 *      Author: Allan
 */

#ifndef EXAMPLE_CO2SOLUBILITY_H_
#define EXAMPLE_CO2SOLUBILITY_H_

void Example_CO2Solubility()
{
	const double Pmin  = 10.0;
	const double Pmax  = 100.0;
	const double dP    = 0.5;
	const double nC    = 3.0; 
	
	GeoManager manager;
	
	manager.SetAqueousSpecies({"H2O", "H[+]", "Na[+]", "Cl[-]", "HCO3[-]", "CO2(aq)", "OH[-]"});
	//manager.SpeciateAqueousSpecies({"H2O", "Na[+]", "Cl[-]", "CO2(g)"});
	manager.SetGaseousSpecies({"H2O(g)", "CO2(g)"});
	
	GeoSystem system = manager.AssembleGeoSystem();
	
	GeoEquilibrator equilibrator(system);
	
	GeoState state(system);
	
	double P = Pmin;

	cout << setprecision(6) << fixed;
	
	cout << "# Column 1:  Pressure [bar]" << endl;
	cout << "# Column 2:  Solubility CO2 [molal] (T =  40 C, mNaCl = 3.997)" << endl;
	cout << "# Column 3:  Solubility CO2 [molal] (T =  60 C, mNaCl = 3.997)" << endl;
	cout << "# Column 4:  Solubility CO2 [molal] (T =  80 C, mNaCl = 4.001)" << endl;
	cout << "# Column 5:  Solubility CO2 [molal] (T = 120 C, mNaCl = 4.003)" << endl;
	cout << "# Column 6:  Solubility CO2 [molal] (T = 140 C, mNaCl = 4.003)" << endl;
	cout << "# Column 7:  Solubility CO2 [molal] (T = 160 C, mNaCl = 4.003)" << endl;
	cout << "# Column 8:  Solubility CO2 [molal] (T =  40 C, mNaCl = 5.999)" << endl;
	cout << "# Column 9:  Solubility CO2 [molal] (T =  60 C, mNaCl = 5.999)" << endl;
	cout << "# Column 10: Solubility CO2 [molal] (T =  80 C, mNaCl = 5.999)" << endl;
	cout << "# Column 11: Solubility CO2 [molal] (T = 120 C, mNaCl = 5.999)" << endl;
	cout << "# Column 12: Solubility CO2 [molal] (T = 140 C, mNaCl = 5.999)" << endl;
	cout << "# Column 13: Solubility CO2 [molal] (T = 160 C, mNaCl = 5.999)" << endl;
	cout << endl;
	
	while(P <= Pmax)
	{
		state.SetPressure(P);
			cout << setw(15) << left << P;
		
		state.SetTemperature(40.0);
			equilibrator.ResetConditions();
			equilibrator.ImposeChargeBalance();
			equilibrator.ImposeSpeciesMass("H2O", 55.508);
			equilibrator.ImposeElementMass("Na", 3.997);
			equilibrator.ImposeElementMass("Cl", 3.997);
			equilibrator.ImposeElementMass("C", nC);
			equilibrator.Equilibrate(state);
			cout << setw(15) << left << nC - state.Get("Moles", "CO2(g)");
			
		state.SetTemperature(60.0);
			equilibrator.ResetConditions();
			equilibrator.ImposeChargeBalance();
			equilibrator.ImposeSpeciesMass("H2O", 55.508);
			equilibrator.ImposeElementMass("Na", 3.997);
			equilibrator.ImposeElementMass("Cl", 3.997);
			equilibrator.ImposeElementMass("C", nC);
			equilibrator.Equilibrate(state);
			cout << setw(15) << left << nC - state.Get("Moles", "CO2(g)");
			
		state.SetTemperature(80.0);
			equilibrator.ResetConditions();
			equilibrator.ImposeChargeBalance();
			equilibrator.ImposeSpeciesMass("H2O", 55.508);
			equilibrator.ImposeElementMass("Na", 4.001);
			equilibrator.ImposeElementMass("Cl", 4.001);
			equilibrator.ImposeElementMass("C", nC);
			equilibrator.Equilibrate(state);
			cout << setw(15) << left << nC - state.Get("Moles", "CO2(g)");
			
		state.SetTemperature(120.0);
			equilibrator.ResetConditions();
			equilibrator.ImposeChargeBalance();
			equilibrator.ImposeSpeciesMass("H2O", 55.508);
			equilibrator.ImposeElementMass("Na", 4.003);
			equilibrator.ImposeElementMass("Cl", 4.003);
			equilibrator.ImposeElementMass("C", nC);
			equilibrator.Equilibrate(state);
			cout << setw(15) << left << nC - state.Get("Moles", "CO2(g)");
			
		state.SetTemperature(140.0);
			equilibrator.ResetConditions();
			equilibrator.ImposeChargeBalance();
			equilibrator.ImposeSpeciesMass("H2O", 55.508);
			equilibrator.ImposeElementMass("Na", 4.003);
			equilibrator.ImposeElementMass("Cl", 4.003);
			equilibrator.ImposeElementMass("C", nC);
			equilibrator.Equilibrate(state);
			cout << setw(15) << left << nC - state.Get("Moles", "CO2(g)");
			
		state.SetTemperature(160.0);
			equilibrator.ResetConditions();
			equilibrator.ImposeChargeBalance();
			equilibrator.ImposeSpeciesMass("H2O", 55.508);
			equilibrator.ImposeElementMass("Na", 4.003);
			equilibrator.ImposeElementMass("Cl", 4.003);
			equilibrator.ImposeElementMass("C", nC);
			equilibrator.Equilibrate(state);
			cout << setw(15) << left << nC - state.Get("Moles", "CO2(g)");
		
		state.SetTemperature(40.0);
			equilibrator.ResetConditions();
			equilibrator.ImposeChargeBalance();
			equilibrator.ImposeSpeciesMass("H2O", 55.508);
			equilibrator.ImposeElementMass("Na", 5.999);
			equilibrator.ImposeElementMass("Cl", 5.999);
			equilibrator.ImposeElementMass("C", nC);
			equilibrator.Equilibrate(state);
			cout << setw(15) << left << nC - state.Get("Moles", "CO2(g)");
			
		state.SetTemperature(60.0);
			equilibrator.ResetConditions();
			equilibrator.ImposeChargeBalance();
			equilibrator.ImposeSpeciesMass("H2O", 55.508);
			equilibrator.ImposeElementMass("Na", 5.999);
			equilibrator.ImposeElementMass("Cl", 5.999);
			equilibrator.ImposeElementMass("C", nC);
			equilibrator.Equilibrate(state);
			cout << setw(15) << left << nC - state.Get("Moles", "CO2(g)");
			
		state.SetTemperature(80.0);
			equilibrator.ResetConditions();
			equilibrator.ImposeChargeBalance();
			equilibrator.ImposeSpeciesMass("H2O", 55.508);
			equilibrator.ImposeElementMass("Na", 5.999);
			equilibrator.ImposeElementMass("Cl", 5.999);
			equilibrator.ImposeElementMass("C", nC);
			equilibrator.Equilibrate(state);
			cout << setw(15) << left << nC - state.Get("Moles", "CO2(g)");
			
		state.SetTemperature(120.0);
			equilibrator.ResetConditions();
			equilibrator.ImposeChargeBalance();
			equilibrator.ImposeSpeciesMass("H2O", 55.508);
			equilibrator.ImposeElementMass("Na", 5.999);
			equilibrator.ImposeElementMass("Cl", 5.999);
			equilibrator.ImposeElementMass("C", nC);
			equilibrator.Equilibrate(state);
			cout << setw(15) << left << nC - state.Get("Moles", "CO2(g)");
			
		state.SetTemperature(140.0);
			equilibrator.ResetConditions();
			equilibrator.ImposeChargeBalance();
			equilibrator.ImposeSpeciesMass("H2O", 55.508);
			equilibrator.ImposeElementMass("Na", 5.999);
			equilibrator.ImposeElementMass("Cl", 5.999);
			equilibrator.ImposeElementMass("C", nC);
			equilibrator.Equilibrate(state);
			cout << setw(15) << left << nC - state.Get("Moles", "CO2(g)");
			
		state.SetTemperature(160.0);
			equilibrator.ResetConditions();
			equilibrator.ImposeChargeBalance();
			equilibrator.ImposeSpeciesMass("H2O", 55.508);
			equilibrator.ImposeElementMass("Na", 5.999);
			equilibrator.ImposeElementMass("Cl", 5.999);
			equilibrator.ImposeElementMass("C", nC);
			equilibrator.Equilibrate(state);
			cout << setw(15) << left << nC - state.Get("Moles", "CO2(g)");
			
		cout << endl;
		
		P += dP;
	}
}

#endif /* EXAMPLE_CO2SOLUBILITY_H_ */
